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We investigate extinction dynamics in the paradigmatic model of two competing species A and B 
that reproduce (A — > 2A, B 2B), self-regulate by annihilation (2A — > 0, 2B — > 0), and compete 
(A + B — > A, A + B — > B). For a finite system that is in the well-mixed limit, a quasi-stationary 
state arises which describes coexistence of the two species. Because of discrete noise, both species 
eventually become extinct in time that is exponentially long in the quasi-stationary population size. 
For a sizable range of asymmetries in the growth and competition rates, the paradoxical situation 
arises in which the numerically disadvantaged species according to the deterministic rate equations 
survives much longer. 
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In the paradigmatic two-species competition model, a 
population is comprised of distinct species A and B, each 
of which reproduce and self regulate by intraspecies com- 
petitive reactions. In addition, interspecies competitive 
reactions occur, which are deleterious to both species 
For large, well-mixed populations, the dynamics can be 
accurately described by deterministic rate equations. For 
finite systems, however, fluctuations in the numbers of in- 
dividuals ultimately lead to extinction, in stark contrast 
to the rate equation predictions. 

In this work, we investigate how asymmetric inter- 
species competition influences the extinction probability 
of each species. In a finite ecosystem, extinction arises 
naturally when multiple species compete for the same re- 
sources. In such an environment, one species often dom- 
inates, while the others become extinct HHil, a feature 
that embodies the competitive exclusion principle. A re- 
lated paradigm appears in the context of competing par- 
asite strains that exploit the same host population, or in 
the fixation of a new mutant allele in a haploid popula- 
tion whose size is not fixed Q- 

With asymmetric interspecies competition, we uncover 
the surprising feature that deterministic and stochastic 
effects, which originate from the same elemental reac- 
tions, act oppositely. For sizable asymmetry ranges in the 
growth and competition rates, the situation arises where 
the combined effects of the elemental reactions leads to 
one species being numerically disadvantaged at the mean- 
field level, despite its interspecies competitive advantage, 
but this competitive advantage dominates the other re- 
action processes at the level of large deviations. Thus 
the outcompeted and less abundant species has a higher 
long-term survival probability: "survival of the scarcer" . 

Model: Asymmetric competition of two species A and B 
is defined by the reactions: 
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The first line accounts for reproduction, the second for in- 
traspecies competition, and the last for interspecies com- 
petition. Here K is the environmental carrying capacity, 
which sets the size of the overall population, e quantifies 
the severity of the competition, while g and a quantify 
the asymmetries in the growth and interspecies compe- 
tition rates, respectively. In our presentation, we focus 
on the limit K ^> 1. While a general model should also 
contain asymmetry in the intraspecies competition rate, 
no new phenomena arise by this generalization; for sim- 
plicity, we study the model defined by Eqs. (QJ. 

To probe extinction in two-species competition, we fo- 
cus on P m , n {t), the probability that the population con- 
sists of m > As and n > Bs at time t. In the limit 
of a perfectly-mixed population, the stochastic reaction 
processes in ([I]) lead to P m , n (t) evolving by the master 
equation 
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Here E and F are the raising and lowering operators [8( 
for species A and B, respectively; viz. WP mtn = P m +i,n 
and P m n P ln n +j. 

Deterministic Rate Equations: First we focus on the av- 
erage population sizes (m) = J2 m n m Pm,n and (n) = 
J2 m n n Pm,n- From ([2]), the evolution of these quantities 
is given by 
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Here we neglect correction terms of the order of 1/K 
and, more importantly, neglect correlations by assuming 
that (to 2 ) = (m) 2 , (n 2 ) = (n) 2 , and {ran) = (m)(n). 
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We restrict ourselves to the parameter range ae < g < 
1/e, which guarantees that the fixed point corresponding 
to coexistence of both species is stable. The four fixed 
points of the rate equations © are then: 
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unstable node, 
saddles, (4) 

stable node. 



If the initial populations of both species are non-zero, 
they are quickly driven to the stable node (Fig. [1} that 
describes the steady-state populations in the mean-field 
limit. The relaxation time toward the stable node, r r , is 
independent of K. These steady-state populations of the 
two species are equal when 
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For g < g* , the i?-population is scarcer. Naively, the 
scarcer population should be more likely become extinct 
first. However, as we shall show, a proper account of the 
fluctuations that stem from the underlying elemental re- 
actions themselves leads to a radically different outcome. 




FIG. 1. Schematic flow diagram in asymmetric two-species 
competition for weak competition in the mean field. The un- 
stable node, the saddles, and the stable node are shown as 
open, hatched, and solid, respectively. 



Extinction: The mean-field picture is incomplete because 
fluctuations of the population sizes about their fixed- 
point values are ignored. For large populations (corre- 
sponding to carrying capacity K 3> 1), these fluctuations 
are typically small. Thus the populations achieve a quasi- 
stationary state where the two species coexist. This 
state is stable in the mean-field description (heavy dot 
in Fig. [T]). However, an unlikely sequence of deleterious 
events eventually occurs that ultimately leads one popu- 
lation, and then the other, to extinction. After the first 
extinction, the remaining population settles into another 
quasi-stationary state around one of the single-species 
fixed points (m*,n*) = (K, 0) or (0, Kg). Eventually a 
large fluctuation drives the remaining species to extinc- 
tion. This second extinction time is typically much longer 
[by a factor that scales as exp(const. x K)] than the first 
time, because the remaining species does not suffer inter- 
species competition. Once a species is extinct, there is 



no possibility of recovery since there is no replenishment 
mechanism. 

The question that we address is: which species typi- 
cally goes extinct first? The answer is encoded in the 
dynamics of the two-species probability P m ,n(t). During 
the initial relaxation stage, a quasi-stationary probability 
distribution is quickly reached (Fig. [2]). The probability 
distribution is sharply peaked at the stable fixed point 
of the mean-field theory This probability slowly "leaks" 
into localized regions near each of the single-species fixed 
points (m*,n*) = (K,0) and (0, Kg). Thus two sharply- 
peaked single-species distributions start to form. If the 
(K, 0) peak grows faster then the B species is more likely 
to go extinct first. Similarly, a faster growing (0, Kg) 
peak means A extinction is more likely. Eventually, the 
probability distribution that is localized at one of the two 
single-species fixed points slowly leaks toward the fixed 
point (0,0) that corresponds to complete extinction Q. 

To determine extinction rates, it is helpful to define 
Pa, Pb, and V$ as the respective probabilities that 
species A is extinct, species B is extinct, or neither is 
extinct at time t [lOj. (Being interested in times much 
shorter than the expected extinction time of both species, 
we can neglect the probability of the latter process.) By 
definition, these extinction probabilities are 
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these satisfy Pa + Pb + P<j> = lj up to an exponen- 
tially small correction that stems from the process where 
both species become extinct simultaneously. In the limit 
K 1 and for times much greater than the relaxation 
time scale r r , the sums in Eqs. (j6|) are dominated by 
contributions from values of to and n that are close to 
the single-species and coexistence fixed points. Moreover, 
these extinction probabilities evolve according to a set of 
effective coupled equations 



Pa = RaP^ , 
P B = RbP^ , 
P^ = -{R a + Rb)P^- 



(7) 



that define Ra and Rb as the respective extinction rates 
for species A and species B. Solving these equations yields 
the time dependence of the extinction probabilities 
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with P. = Ra + Rb- To determine Ra and Rb, we follow 
the evolution of the eigenstate of the master equation 
([2]) that determines the leakage of probability from the 
vicinity of the coexistence point: 
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and 1Z is the third-lowest positive non-trivial eigenvalue 
of the operator H . The two still-smaller positive non- 
trivial eigenvalues correspond to the much slower decay 
of quasi-stationary single-species states and play no role 
in the dynamics of the first extinction event. There is also 
a trivial eigenvalue that corresponds to the final state of 
complete extinction. 

Combining Eq. ((2]) with ©-((9]), we obtain the follow- 
ing expression for the extinction rate of the A species: 



Ra = j£ (e n U.i >n + n 2 ,„) 

n>0 
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As expected, the extinction rate for As involves two pro- 
cesses: (i) elimination of the last remaining A via compe- 
tition with Bs and (ii) annihilation of the last remaining 
pair of As. Similarly, 
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To calculate Ra and Rb, we therefore need to evaluate 
the small-population-size tails of H m „. This task can 
be achieved by applying a variant of Wentzel-Kramers- 
Brillouin (WKB) approximation, that was pioneered in 
Refs. [llMl4| . and was applied more recently to pop- 
ulation extinction, in particular, for stochastic two- 
population systems [l(| IB"22|. The WKB ansatz for 
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where x = m/K and y = n/K are treated as continuous 
variables. Substituting Eq. (|12[) into (|llap and assuming 
K 3> 1, gives, to lowest order in 1/K 



R A ~e~ KSA , R B ~e- K ° B , (13) 

where S A = S(0,g) and S B = 5(1,0), (see Refs. 
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17[). Thus as K ^> 1, the eventual extinction probabili- 
ties in ([5]) simply become (up to pre-exponential factors 
that depend on K), as 



V A (t = oo) = l-V B (t = oo) 
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To determine the extinction probabilities explicitly, we 
therefore need Sa and Sb- To this end, we substitute the 
WKB ansatz ([12]) in Eq. ([10]) and Taylor expand S{x,y) 
to lowest order in 1/K. After some algebra, we obtain an 
effective Hamilton- Jacobi equation H(x,y,d x S,d y S) = 
—TZ, with the Hamiltonian 



H(x,y, Px , Py ) = x(e p *-l) + gy(e p y -1) 



.fL( e -ap«_i) 



(15) 
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FIG. 2. Quasi-stationary probability distributions for species 
A, P m = z^n Prn,n, and species B, P„ = J2 m p m,n- Parame- 
ters are K = 100, e = 0.9, g = 0.45, and a = 0. Symbols are 
simulation results, while the solid curve is WKB approxima- 
tion for the B species distribution. 



Here p x = d x S and p y = d y S are the canonical momenta 
that are conjugate to the "coordinates" x and y. Corre- 
spondingly, S(x, y) is the classical action of the system. 

Since we expect that —1Z (which now has the meaning 
of energy in this Hamiltonian system) is expected to be 
exponentially small, we set it to zero. The Hamiltonian 
equations of motion x = dH/dp x , p x — —dH/dx, etc., 
have six finite zero-energy fixed points, and three more 
fixed points where one or both momenta are at minus 
infinity. Only three of the fixed points, however, turn 
out to be relevant for answering our question about which 
species typically goes extinct first. These are: 
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F A = (0,.g,ln( 5 e),0) , 
F B = (l,0,0,ln(ae/ 3 )) 
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A straightforward way to determine S A and Sb is by 
calculating the action along the activation trajectories. 
These are zero-energy, but non-zero-momentum trajec- 
tories of the Hamiltonian system ([15]) that go from F$ to 
F A , and from Fa to Fb, respectively. These actions are 



S A = 



{p x dx + Pydy) , 



(17) 



and similarly for Sb- In general, these activation 
trajectories — separatrices, or instantons — cannot be cal- 
culated analytically because of the lack of an integral of 
motion that is independent of the energy. However, for 
e <C 1 a pcrturbative solution for these trajectories is 
possible. 

As a preliminary, we outline how to calculate the ac- 
tion for the special case e = 0, which corresponds to 
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two uncoupled species. Here the zero-energy activation 
trajectories can be easily found. For the — > Fa sep- 
aratrix, the B species is unaffected by A extinction so 
(y,p y ), which correspond to the coordinates of the Bs, 
remains constant throughout the evolution. As a result, 
a parametric form of the F^ — > Fa separatrix is 
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with y = g and p y = throughout. Similarly, for the 
F$ — > Fb separatrix one obtains 
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with x — 1 and p x — throughout. Substituting the 
trajectories given in (|18[) into (|17[) and performing the 
integration by parts gives Sa = 2(1 — In 2) and Sb = 
2g(l -ln2) $i\24\. 

For weak interspecies competition (e <C 1), we can 
calculate the corrections to the actions to first order 
in e. For this purpose, we split the Hamiltonian f| 15[) 



into unperturbed and perturbed parts, H = Hq + eH\, 
and similarly expand the action as S = So + eSi + . . . . 
Following [HI [3, HH, the correction to the action is 

Sl = I-oo Hl \- x ( t ^y( t ^P x ( t ^Pv( t ^ dt ^ wnere the integral 
is evaluated along the unperturbed trajectories given by 
Eqs. (fT8|) . Performing this integral for Si yields the cor- 
rected actions 
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2(1 - In 2) -e(2c/ln2), 
2 5 (1 - In 2) -e(2aln2) 
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Equations ([131 and (fT9|) give the analytic expression 
for the extinction probabilities of each species for weak 
interspecies competition. Using Eq. (|19[) and imposing 
the condition Sa = Sb from Eq. (fl4)) . we obtain the 
following condition for equal extinction probability for 
both species: 
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where c = [(In 2) _1 - 1] ~ . The predictions of Eqs. (TH)) 
and (Tig)) are in good agreement with our simulation re- 
sults (Fig.©. 

Phase Diagram: Comparing Eqs. ([5|) and (|20|) . one sees 
that there is a sizable region in the a-g parameter space 
where one species has a smaller quasi-stationary popu- 
lation and yet an (exponentially) smaller probability to 
first become extinct. As an illustration, Fig. [3] shows the 
probability for B to become extinct first for fixed a and 
e. We also produced analogous curves as Fig. [3] at many 
values of a. From the value of g at which the extinction 
probabilities are equal, we infer the phase diagram shown 
in Fig. 01 Simulations at larger values of e yield the same 
qualitative phase diagram. 




FIG. 3. Probability that species B first becomes extinct as 
a function of growth rate asymmetry g for a — 0, e = 0.1, 
K = 40. The curve is the prediction of Eq. (|14|l . with actions 
given by (| ltJf) . while circles are simulation results. In the 
hatched region, the quasi-stationary population of B is less 
than that of A, but Bs are less likely to become extinct first. 




FIG. 4. Phase diagram for e = 0.1 showing loci of equal 
quasi-stationary population sizes (dashed) and equal extinc- 
tion probabilities from Eq. (|20[) (solid). Circles indicate sim- 
ulation results. In the hatched region, As are more numerous 
in the quasi-stationary state but are more likely to become ex- 
tinct first. In the cross-hatched region Bs are more numerous 
but are more likely to become extinct first. 



Conclusion: In two-species competition, interspecies 
competitive asymmetry leads to the unexpected phe- 
nomenon of "survival of the scarcer" . The very same 
elemental reactions that lead to a disadvantage in the 
quasi-stationary population size of one species within a 
deterministic mean-field theory, may also give this species 
a great advantage in its long-term survival when fluctu- 
ation effects are properly accounted for. 
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